void landscapediff(double *h,double *old_h,double k,double dt) {
    double D;
    double *a,*b,*c,*d,*bdot,*ddot;
    
    a = calloc(Sbase,sizeof(double));
    b = calloc(Sbase,sizeof(double));
    c = calloc(Sbase,sizeof(double));
    d = calloc(Sbase,sizeof(double));
    bdot = calloc(Sbase,sizeof(double));
    ddot = calloc(Sbase,sizeof(double));
    
    for (s=1;s<Sbase-1;s++) {
        D=dx[s]+dx[s-1];
        a[s]=2*k/(D*dx[s-1]);
        b[s]=-(2*k/(D*dx[s-1])+2*k/(D*dx[s])+1/dt);
        c[s]=2*k/(D*dx[s]);
        d[s]=-old_h[s]/dt;
    }
    d[0]=0; d[Sbase-1]=0;
    b[0]=1; c[0]=-1;
    b[Sbase-1]=1; a[Sbase-1]=-1;
    tridiag(h,a,b,c,d,bdot,ddot,Sbase);
    
    free(a);
    free(b);
    free(c);
    free(d);
    free(bdot);
    free(ddot);
}


//as above, but variable diffusivity
void landscapediff1(double *h,double *old_h,double *k,double dt) {
    double alphaplus,alphaminus;
    double *a,*b,*c,*d,*bdot,*ddot;
    
    a = calloc(Sbase,sizeof(double));
    b = calloc(Sbase,sizeof(double));
    c = calloc(Sbase,sizeof(double));
    d = calloc(Sbase,sizeof(double));
    bdot = calloc(Sbase,sizeof(double));
    ddot = calloc(Sbase,sizeof(double));
    
    for (s=1;s<Sbase-1;s++) {
        alphaplus=1/(dx[s]*(dx[s-1]+dx[s]))*(k[s]+k[s+1]);
        alphaminus=1/(dx[s-1]*(dx[s-1]+dx[s]))*(k[s]+k[s-1]);
        a[s]=alphaminus;
        b[s]=-(alphaminus+alphaplus+1/dt);
        c[s]=alphaplus;
        d[s]=-old_h[s]/dt;
    }
    
    d[0]=0; d[Sbase-1]=0;
    b[0]=1; c[0]=-1;
    b[Sbase-1]=1; a[Sbase-1]=-1;
    tridiag(h,a,b,c,d,bdot,ddot,Sbase);
    
    free(a);
    free(b);
    free(c);
    free(d);
    free(bdot);
    free(ddot);
}



void surfaceupdate(int newsurface) {
    double *newtopo,*oldtopo,*oldsedipool,*k_var;
    double k;
    double *oldsurfaces;
    double topo;
    double slope, interc;
    double sedipoolincrement;
        
    if (newsurface) { //create a new surface, allocate new memory and copy old surfaces
        oldsurfaces=surfaces;
        surfaces=calloc((Nsurfaces+1)*Sbase,sizeof(double));
        for (e=0;e<Nsurfaces*Sbase;e++) surfaces[e]=oldsurfaces[e];
        Nsurfaces++;
        free(oldsurfaces);
        for (s=0;s<Sbase;s++) surfaces[INsur(s,Nsurfaces-1)]=surfaces[INsur(s,Nsurfaces-2)]; //copying the old top surfaces to the new top surface
    }
    
    newtopo=calloc(Sbase, sizeof(double));
    oldtopo=calloc(Sbase, sizeof(double));
    oldsedipool=calloc(Sbase, sizeof(double));
    k_var=calloc(Sbase, sizeof(double));
    
    sealevel=0.5*(surfaces[INsur(0,3)]+surfaces[INsur(Sbase-1,3)]);
    
    for (s=0;s<Sbase;s++) sedipool[s]=1e4;



  
    for (s=0;s<Sbase;s++) {
        if (surfaces[INsur(s,Nsurfaces-1)]>sealevel) {  //sedimentation criterium
            surfaces[INsur(s,Nsurfaces-1)]-=sedipool[s];
            sedipool[s]=0; //emptying sedipool
            if (surfaces[INsur(s,Nsurfaces-1)]<sealevel) { //if overflow - return surplus to sedipool
                sedipool[s]=sealevel-surfaces[INsur(s,Nsurfaces-1)]; surfaces[INsur(s,Nsurfaces-1)]=sealevel; }}}
    
    
/*
    for (s=0;s<Sbase;s++) {
        surfaces[INsur(s,Nsurfaces-1)]=newtopo[s];
        for (sur=Nsurfaces-2;sur>=0;sur--) {
            if (surfaces[INsur(s,Nsurfaces-1)]>surfaces[INsur(s,sur)]) surfaces[INsur(s,sur)]=surfaces[INsur(s,Nsurfaces-1)];
                else break; }
    }
    */
    
    //updating markers according to the new surfaces
    for (e=0;e<N;e++) {
        s=mcs[e]; t=mct[e];
        slope=(surfaces[INsur(s+1,Nsurfaces-1)]-surfaces[INsur(s,Nsurfaces-1)])/(x[s+1]-x[s]);
        interc=surfaces[INsur(s,Nsurfaces-1)]-slope*x[s];
        topo=slope*mx[e]+interc; //local topography
        if (topo<mz[e]) { //possible sedimentation
            if (mtype[e]==0) { //check for air to sediment conversion
                mtype[e]=5; //WARNING: seditype may vary
                mstrainxx[e]=1e-20; mstrainxz[e]=1e-20; //decrease strain for the next time step
            }
        }
        else { //possible erosion
            if (mtype[e]>0) {
                //mtype[e]=0; //mstressxx[e]=0; mstressxz[e]=0;
            }
        }
    }
    free(newtopo);
    free(oldtopo);
    free(oldsedipool);
    free(k_var);
}
